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Abstract. 

We study how the Newton-GMRES iteration can enable dynamic simulators (time-steppers) to 
perform fixed-point and path-following computations. For a class of dissipative problems, whose 
dynamics are characterized by a slow manifold, the Jacobian matrices in such computations are 
compact perturbations of the identity. We examine the number of GMRES iterations required for 
each nonlinear iteration as a function of the dimension of the slow subspace and the time-stepper 
reporting horizon. In a path-following computation, only a small number (one or two) of additional 
GMRES iterations is required. 
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1. Introduction. In studying the dynamic behavior of evolution equations, 
(1.1) du/dt = f(u;X), 

a computational modeler typically chooses between two paths: the first is developing 
a dynamic simulator for the process; the second is developing algorithms to locate 
particular features of the long-term dynamics of the process, such as steady states 
or limit cycles. The first path typically gives rise to initial value problems, and the 
corresponding codes are dynamic simulators - we will call them "time-steppers" since, 
given the state of the system at a moment in time they produce (an approximation 
of) the state of the system at a later moment. The second path typically gives rise to 
fixed point algorithms for solving coupled nonlinear algebraic equations; these may 
be the steady state equations themselves, or an augmented set arising in a continua- 
tion/bifurcation context. The Recursive Projection Method of Shroff and Keller [24] 
(see also the Adaptive Condensation of Jarausch and Mackens [12, 13] for symmetric 
problems and the Newton-Picard work of Lust et al. [20] ) is an example of an algo- 
rithm that, in some sense, connects the two paths. The main idea is to construct a 
computational superstructure that designs and combines several calls to an existing 
("legacy") time-stepper, effectively turning it into a fixed point solver for 

u - $t(w; A) = 0, 

where <£>t is the result of integration of (|l.lfl with initial condition u for time T (see 
the review of Tuckerman and Barkley [30]). Here the action of the linearization of the 
time-stepper is estimated in a matrix-free fashion by the integration of appropriately 
chosen nearby initial conditions. 

Over the last few years, this matrix-free computational enabling technology has 
found new applications in the field of multiscale computations. In current modelling 
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practice, dynamic models are often constructed at a microscopic/stochastic level of de- 
scription (e.g. molecular dynamics, kinetic Monte Carlo or Lattice Boltzmann codes); 
the closure required to obtain explicit, macroscopic "effective" equations is not avail- 
able in closed form. Algorithms like RPM can then be "wrapped around" appropri- 
ately initialized ensembles of short bursts of microscopic simulations, and solve the 
effective equation without ever obtaining it in closed form. Several examples of such 
"coarse equation-free" computation have now been explored, and the mathematical 
underpinnings of the approach are being extensively studied [9,19,26,28]. 

In this paper we study the Newton-GMRES iteration as a computational "wrap- 
per" around a legacy time-stepper. This wrapper enables the computation and contin- 
uation of fixed points of the time-T map of the time-stepper (i.e. steady states of the 
corresponding dynamical equations). It is useful, for purposes of discussion, to con- 
sider that the time-stepper is available as an input-output black box (an executable) 
which cannot be modified. 

The paper is organized as follows: We first review certain properties of GMRES 
in the context of problems whose linearization is a compact perturbation of the iden- 
tity. We argue that time-steppers for a class of dissipative problems, whose long-term 
dynamics lie on a low-dimensional, slow manifold, may fit this description; we then 
proceed to examine Newton-GMRES convergence and number of iterations for fixed 
point computation of such time-steppers in a continuation context. Our numerical ex- 
amples are the Chandrasekhar H-equation [5] and the time-stepper of a discretization 
of a reaction-diffusion problem known to possess a low-dimensional inertial mani- 
fold [8,27]. We illustrate the effect, on the GMRES iterations, of the time-stepper 
reporting horizon both for fixed point and for pseudo-arclength continuation compu- 
tations. We conclude with a brief discussion and thoughts about extensions of this 
approach. 

2. Convergence Analysis. 

2.1. GMRES Preliminaries. We will use Newton-GMRES to solve the non- 
linear fixed point equation 

(2.1) u-<f> T (u;\) = F(u) = 0. 

In what follows we will work in R N and use the usual Euclidean norm; the idea is that 
equation (|2.1|l arises from integrating a set of ODEs, possibly from the discretization 
of a partial differential equation on a given mesh. We will also assume (and the 
consequences of this will become apparent below) that the long-term dynamics of the 
discretized PDE occur on an attracting "slow manifold" of low dimension p << N; 
p will be assumed independent of mesh refinement (i.e. it will remain constant as N 
increases.) 

GMRES is an iterative method for solving linear systems Ax = b in R N . The 
kth GMRES iteration minimizes the residual r over xo + fCk, where xo is the initial 
iterate and ICk is the the kth Krylov subspace 

JCk = span{r , Ar , . . . , A ~Vo}. 

A consequence [15,23] of the minimization is that 

(2.2) ||r fe || = min ||p(A)r || 

where Vk is the space of kth degree residual polynomials, i.e. polynomials of degree 
k such that p(Q) = 1. 
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We will apply l|2.2[l to a special class of problems where 
(2.3) A = I - K + E. 

In JOJ 

• I — K is nonsingular, 

• K = P]jKP]j, where Pjj is an orthogonal projection onto a space D of 
dimension p << N, and 

• E is a matrix with small norm. 

We will analyze the performance of GMRES in a way different from the eigenvalue- 
based approach used for diagonalizable matrices [10, 15, 18, 29]. For the class of prob- 
lems of interest here, we can prove a convergence result directly using methods similar 
to those in [3,4, 7, 17]. The result in this paper has stronger and sharper convergence 
rates. This result carries through in the infinite-dimensional case also, using the L 2 
norm for the corresponding function spaces. 

Theorem 2.1. Let A be given by (|2.3|l . Then there exists C such that for all 
m > 1, 

(2-4) ||r m(p+1) ||<C m ||i;|r. 

Proof. Let pc be the characteristic polynomial of / — K. Since D has dimension 
p, pc has degree p + 1 . Clearly 

p c (A) =p c (I-K) + A = A. 

by the Cayley-Hamilton theorem. Moreover, there is C > such that 

l|A||<C||i5||. 

In fact, if 

P+i 

p{z) = l + ^ 7fc z fc 
fc=i 

then we may use 

p+i 

C = Y, k \lk\\\I - P D KP D \\ k - x + 0(||i?|| 2 ). 

fc=i 

Hence 

\\P(A)\\<C\\E\\. 

This is (El for m = 1. 
Define 

p{z) =p c (z)/pc(0). 
Clearly p <E V p +i- Hence, by (|2~!>|) 

(2.5) \\r m{p+1) \\<\\p m (A)r \\<C m \\E\n\r \\, 
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as asserted. □ 

The estimate in <|2.5(1 does not depend on the eigenvalues of A, nor is there any 
requirement that A be normal or even diagonalizable. 

As is standard, a GMRES iteration is terminated when 

Nil < v\\ro\\, 

where r\ is an user-defined parameter. In the context of this paper \\E\\ is well below 
the termination tolerance r, so we conclude that the iteration will terminate in at 
most p + 1 iterations. In the general case, of course, more cycles could be required. 

2.2. Time-Steppers and Steady State Solutions. Let $(it) denote the out- 
put of the time-stepper with time step T and initial data u; we have dropped the 
subscript T of <E> for convenience. We seek to solve 

(2.6) F(u) = u - = 0, 

to find a steady state solution. Consider the structure of the eigenvalues fii = exp(tTiT) 
of given the structure of the eigenvalues <7i of the linearization of the original 
problem. We assume that there exists a significant gap between a few Ui close to zero 
(p of them, to be exact) and the remaining large negative Uj. The p small norm Oi do 
not have to be stable - they can also be slightly unstable, so that the corresponding 
Hi are close to the unit circle. 

We then assume that T is large enough so that there is a space D of low dimension 
p, so that 

$(u) = P D <f>(P D u) + E(u) 

and E and its Jacobian E u are small. For a problem where some of the p small Oi 
are unstable, this large time T should not be so large that the unstable modes will 
cause the trajectory to move significantly away from the fixed point. In such cases 
the proper choice of T is a somewhat delicate matter [19,24]. 

The equation for the Newton step from a current iterate u c is 

(2.7) F u (u c )s = s - $ u (u c )s = -F(u c ), 

here F u and <&„ are the Jacobians of F and $. 

A Newton-GMRES [1,2,15] method solves l(2~7|l with GMRES, terminating the 
linear (or inner) iteration when 

\\F u (u c )s\\ <Vc\\F(u c )\\ 

where rj c may be changed as the outer (or nonlinear) iteration progresses [6,15,16]. 
For the application here, we assume that ||J5 u (u)|| is much less than any choice of 77 
we make during the iteration. 

The linear system ()2.7|l fits exactly in to the paradigm of § 12.11 with K = 
Pd&u(Pdu)Pd and E = E u (u). Hence we conclude that a GMRES iteration for 
(12 .7(1 will take at most p + 1 iterations to drive the residual to 0(||i?||). 

2.3. Time- Steppers and Continuation. In a continuation context <& depends 
on a parameter A. As is standard [14] we add an additional arclength parameter s to 
obtain the augmented system 



(2.8) 
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In (|2.8p . u and A are approximations to du/ds and dX/ds, which can be obtained in 
several ways [14]. In the calculations reported in §[3|we used the estimate of the slope 
given by the last two points, (uq, Xq) and (it-i, A_i) computed on the branch. 
G is defined on R N+1 and 



G Ul x(u,X) 



F u F x 
u T X 



We seek to show that G Ui x also fits into our paradigm, with p replaced by at most 
p + 2. Hence, at most two additional linear iterations will be needed for each Newton 
iteration of the augmented system. 
We use (|2.6|) to obtain 



G u ,x = 



I - P D $ u (P D u;X)P D 



Pd$x 
X 



-E u 








Let 



where, 



and 



A 



I 
1 



Pd^uPd Pd^x 
-u T 1 — A 



£ = 







-Ex 




D 




This fits the paradigm of § 12.11 To see that, let 

T> = span 
Clearly the range of JC 

R(JC) = span 



D 




c R 



C V. 



N+l 



To apply the results from the previous section, however, we need K, = P-dJCPd, where 
Pp is the orthogonal projector onto D, This is why the extra dimension (ii, 0) T is 
required. 

To see this, let y = (u, fx) C R N+1 be orthogonal to V. This means that y = (w, 0), 
where u> is orthogonal to both D and it. Clearly 



ICy = 



Pd<S>uPdu 



• T 

-u u> 



so IC(I — Px>) = 0. Since R(fC) C V, we have Pp/C = K,. Summarizing 



K = PvlCPv 
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and 

dim (2?) < p + 2. 

The dimension of T> can be taken p + 1 if u is nearly in the range of Pn , i- e. 
(2-9) 11(7-^11=0(11^11 + 11^11). 

In this case, we can let 

l> = span((?U?)W', 



/c = 

and 

£ = 



) ' V 1 



-(Pdu) t 1 - A 



((/ - P D )u) T 

If F u is well conditioned and |A| and ||-Fa|| are not too large, then l|2.9[l holds. To 
see this we differentiate F = and obtain 

F u u + F x X = 0, 

which implies, if F u is nonsingular, that 

(2.10) u = -F^FxX. 

Our assumptions are that \\E U \\ and \\E\\\ are much smaller than \\F U \\ and \\F\\\. 
Hence, if F u is well conditioned, the Banach lemma implies that 

F- 1 = (I-PoQuPD^ + OiWEvW). 

Moreover, 

F x = -P D <f> x + 0(\\E x \\). 
We incorporate this into (|2.1U|) to obtain 

(2.11) u = (I - Pd^uPd^Pd^xX + O (|A|(||£ ||||Fa|| + \\E X \\)) , 

which implies if |A| and \\F X \\ are 0(1). 

Summarizing, if \\E U \\ and \\E X \\ are much smaller that \\F U \\ and ||-Fa|| respec- 
tively, |A| and are O(l), and if F u is well conditioned, then l|2.9|l holds. Near 
folds and bifurcations, F u is singular, and in those cases we may need to require that 
the dimension of V be p + 2. 

3. Numerical Results. All the computations were done using MATLAB ver- 
sion 6.0 Release 12 on a PC with a Pentium4 2.53GHz CPU. 
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3.1. The H-equation. As a first example, we solved a problem which does not 
arise in a time-stepper context, but for which the Jacobian fits our paradigm [22]. 
The solution and continuation problem for a 100-node midpoint rule discretization of 
the Chandrasekhar H-equation 

l 2^ f~{ IH + H j 

were obtained with a Newton-GMRES solver nsoli from [16]. We set the relative 
and absolute tolerances in the solver to 10 -12 . 

2.8 1 1 1 1 1 



2.6- 




1.4 1 1 1 1 1 

0.85 0.9 0.95 1 1.05 

c 

Fig. 3.1. Bifurcation diagram of the H-equation for c close to 1. The parameter c is equal to 
0.9999179 at the point marked by a circle. 

Figure. 13.11 is a bifurcation diagram with respect to the parameter c, showing a 
turning point at c = 1. For all values of c shown in this continuation, the eigenvalues 
of the linearization of the solution are close to 1; only a single one of them changes 
in marked way, ranging from -0.9 on the upper branch through zero at the turning 
point to 0.5 on the lower branch. 

Figure. l3~2l shows the GMRES performance for a representative value of c (0.99991 
79), close to the turning point, marked by a circle in Figure 13.11 We report on 
computations from the continuation itself, where the initial iterate was the standard 
linear predictor, and from a second stand-alone computation where the initial iterate 
was the solution plus a perturbation function 0.05sin(a;). In Figure l3~21 we plot the 
convergence history for Newton-GMRES in the two cases. The convergence rates and 
the costs of the solves are roughly the same. 

The first column of Table. I5~T1 shows the ten eigenvalues farthest from 1 for the 
linearization of each of the two problems. In the continuation case, there is one more 
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-14 1 1 1 1 1 ' 

5 10 15 20 25 30 

function evaluations 

Fig. 3.2. Convergence plot for the steady state and continuation problem of H-equation 



eigenvalue far away from 1 than in the stand-alone case. In both cases the linearization 
clearly fits the pattern of a compact perturbation of the identity. 



Eigenvalues of steady state problem 


Eigenvalues of continuation problem 


0.0207265 


-5.2022448 


0.9424114 


5.2008250 


0.9900391 


0.9746094 


0.9974043 


0.9828252 


0.9992541 


0.9977360 


0.9998164 


0.9991859 


0.9999612 


0.9998098 


0.9999926 


0.9999577 


0.9999987 


0.9999919 


0.9999998 


0.9999985 



Table 3.1 



Eigenvalues for the linearized steady state and continuation problems of the H-equation at 
c = 0.9999179 



3.2. The Chafee-Infante reaction diffusion problem. We then solved the 
steady state and continuation problems for a discretization of a dissipative reaction- 
diffusion PDE in one dimension, the so-called Chafee-Infante problem, 

u t - \u xx + u 3 - u = 0, x e [0, 7r] 
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with boundary conditions it(0, t) — 0,u(ir, t) = 0. We used 201 finite difference 
discretization points. 




_-| i i i i i i i 

1 2 3 4 5 6 7 

lambda 

Fig. 3.3. Bifurcation diagram for the Chafee- Infante reaction diffusion problem. 

Figure . [3 . 31 shows the bifurcation diagram for this problem for a range of param- 
eter values (0 < A < 7) where up to five different spatially structured steady states 
exist. Our computational tests were performed on the upper stable solution branch 
close to A = 2. In this example we studied both the steady state problem (setting the 
right-hand-side of the finite difference equations equal to zero) and the time-stepper 
formulation (our integration routine was odel5s). The steady state (equations) / fixed 
point (time-stepper) and continuation problems were solved with Newton-GMRES 
solver nsoli for various time reporting horizons. The results for A = 2.1386697 are 
shown in Fig. 13.41 Relative and absolute tolerances were chosen to be 10~ 12 . The ini- 
tial guess for the direct solution was chosen to be the true solution plus a perturbation 
function 0.1sin(a;) . 

One thing should be made clear at this point. Using several time steps of an 
implicit integrator, with the concomitant nonlinear solves, is clearly not an efficient 
way of solving a fixed point problem (a single nonlinear solve). The integrator is used 
here as a "legacy code" , a code that one cannot modify. It is in the context of such 
legacy codes that our approach becomes useful, as well as in the case of multiscale 
computations, where the time evolution is performed by a simulator at a different 
level of description (e.g. Lattice-Boltzmann or kinetic Monte Carlo). 

3.3. Clustering of Eigenvalues. Using the Finite Difference Method, the ODE 
system from discretizing the original PDE has a form of ii{t) = f(u),u £ R n . At a 
steady state u* , the matrix f u has n eigenvalues: cr-j, i = 1, . . . , n. If we solve for u* 
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-14 ' ' ' ' 1 1 1 1 ' 

5 10 15 20 25 30 35 40 45 

function evaluations 

Fig. 3.4. Convergence of time-stepper based fixed point computation for the steady state and 
the continuation of the discretized RD problem. Every circle corresponds to a Newton step. A single 
"function evaluation" consists of integration over one time reporting horizon. 

with a time-stepper, the system to be solved is u(0) — $t(w; A) = 0, which also has n 
eigenvalues: 1 — e <TiT , i = 1, . . . , n. 

We have used the forward in time time-stepper to compute both stable and un- 
stable steady states. We have marked the unstable steady state we computed using a 
reporting horizon of T = 0.1 for A = 4.5710239 using forward in time time-stepping. 
When the parameter A is equal to 2.183867, u* is a stable steady state and all the 
eigenvalues are negative. When we increase T, all the eigenvalues 1 — e CTiT are ap- 
proaching ("clustering at") 1. Clustering of eigenvalues is known to be beneficial for 
GMRES performance. Conversely, when we decrease T, eigenvalues start leaving the 
cluster. This results in additional GMRES iterations. 

To quantify the dependence of the performance of the iteration on T we use the 
number of GMRES iterations at the final step in Figure 13.51 When T is large, we 
consistently see 2 linear iterations. As T is reduced, we see an increase in GMRES 
iterations, as expected. In this particular example, when T = 1.78, the number of 
GMRES iterations increases to 3. 

In at attempt to quantify this further we identify a "cluster" of eigenvalues that 
seems to be correlated with the performance of the GMRES iteration. In Fig. 13.61 
we treat those eigenvalues in the interval [0.84 1] as "in the cluster". The number 
of eigenvalues outside the cluster (or smaller than 0.84) and the number of function 
evaluations needed to finish the last Newton step are compared in Tab. 13.21 below. A 
clear, strong correlation emerges. 

3.4. Conclusion. We have extended and sharpened results for the performance 
of GMRES for discretizations of compact fixed point problems to the special class 
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30 40 
function evaluations 



70 



Fig. 3.5. Convergence of time-stepper based fixed point Newton-GMRES computation for dif- 
ferent time reporting horizons. Every circle corresponds to a Newton step, and m is the number of 
function evaluations for the last Newton . 



T=4 
1=2 

T=1 

T=0.5 

T=0.3 

T=0.1 

T=0.07 

T=0.04 

T=0.02 



- o o o 



o o o 



ooooo o o o o ooood 



O O <D 



O O OOGO 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

0.84 



Fig. 3.6. Twenty smallest eigenvalues of the linearized time-stepper at the steady state for 
different time reporting horizons; the dashed line corresponds to the eigenvalue that first leaves the 
cluster when T = 1.78 . 
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T 


4 


2 


1 


0.5 


0.3 


0.1 


0.07 


0.04 


0.02 


Number of 
eigenvalues outside 
the cluster 








1 


2 


3 


6 


7 


10 


14 


Function evaluations 
needed for clustering 
eigenvalues 


2 


2 


2 


2 


2 


2 


2 


2 


2 


Actual 
total function 
evaluations 


2 


2 


5 


6 


7 


11 


8 


12 


16 



Table 3.2 



of problems that arise in a time-stepper context. The key feature leading to the 
enhancement of GMRES performance is the compactification of the spectrum. In 
previous work, a Green's function based reformulation of the steady state problem 
for elliptic PDEs [7] led to a linearization that was a compact perturbation of the 
identity, and an efficient solution via Newton-GMRES. Time-steppers compactify the 
spectrum of the linearization in a natural way, and their properties can be exploited to 
obtain accurate bounds on the convergence rates of the linear iterations in a Newton- 
GMRES continuation. We show that the additional equation in a pseudo-arclength 
formulation of a parameter-dependent family of nonlinear equations adds at most two 
GMRES iterations when the eigenvalues are well separated, and obtain a bound on 
the convergence in terms of the separation of the spectrum and the dimension of the 
slow subspace (associated with a slow/inertial manifold for the dynamics). 

We reported on numerical results that support the theory, first with an integral 
equation, for which we can numerically demonstrate that all but two of the eigen- 
values of the linearization lie in a tight cluster about 1. The second example was 
a parametric study of the steady state of a discretized parabolic partial differential 
equation implemented through time-steppers. 

The "natural" compactification of the spectrum using time-steppers provides a 
natural connection with the performance of matrix-free iterative methods. This com- 
pactification may prove useful in writing computational wrappers, that will accelerate 
the convergence of legacy dynamic simulators to stationary states. We also expect this 
compactification to assist in writing computational wrappers that will assist dynamic 
simulators at a different level of model description (e.g. kinetic Monte Carlo, Brow- 
nian Dynamics or Molecular Dynamics codes [11,21,25]) to locate stationary states 
and perform continuation/bifurcation analysis of macroscopic system observables in 
the so-called equation-free framework [19]. 
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